#### Fixed ####
## Function to evaluate fixed battery scores for each
## person, after asking each specified question
# questions - which questions to choose from
# n - number of questions to use
# model - fit model by grm() on test dataset
# dataset - test dataset
eval_fixed <- function(questions, n, model, dataset){
  scores <- errors <- matrix(NA, nrow = nrow(dataset), ncol = n)
  
  ## this loop evaluates theta after each asked question
  ## by adding question answers into profile one at a time
  for(i in 1:n){
    answer_profiles <- matrix(NA, nrow = nrow(dataset), ncol = ncol(dataset))
    answer_profiles <- data.frame(answer_profiles)
      
    colnames(answer_profiles) <- colnames(dataset)
    answer_profiles[ , questions[1:i]] <- dataset[ , questions[1:i]]

    current_scores <- factor.scores(object = model,
                                    resp.patterns = answer_profiles[ , names(model$X)],
                                    method = "EAP")
    scores[ , i] <- (current_scores$score.dat$z1)
    errors[ , i] <- (current_scores$score.dat$se.z1)
  }
  return(list(scores = scores, errors = errors))
}

#### Adaptive ####
## Function to determine the questions to ask to one
## person given their answers, item parameters, and selection
## and estimation criteria
# row - used for laply indexing
# cat - Cat object used in next item routine
# n - length of battery
# dataset - test dataset
adaptive_one <- function(row, cat, n, dataset){
  cat@answers <- rep(NA, length(cat@guessing))
  answer_progression <- matrix(NA, nrow = n, ncol = ncol(dataset))
  for(i in 1:n){
    next_item <- selectItem(cat)$next.item
    cat@answers[next_item] <- dataset[row, next_item]
    answer_progression[i, ] <- cat@answers
  }
  colnames(answer_progression) <- c(names(cat@difficulty))
  print(answer_progression)
  return(answer_progression)
}

## Function to simulate the question and response profile
## for each respondent in dataset, given specified battery length.
## Returns an array where  each slice is a matrix with rows of 
## respondents columns of questions.  Each slice represents the
## answer profiles after each new question has been asked.
# model - fit model by grm() on training dataset
# n - length of battery
# dataset - test dataset
adaptive_all <- function(model, n, dataset){
  
  cat = grmCat(model, 21) #default grm GHk = 21
  cat@estimation = "EAP"
  cat@selection = "EPV"
  cat@poly = TRUE

  answer_array <- laply(.data = 1:nrow(dataset),
                        .fun = adaptive_one, cat = cat, n = n, dataset = dataset,
                        .parallel=TRUE,
                        .paropts=list(.errorhandling="pass"))

  answer_array <- aperm(answer_array, perm=c(1,3,2))
  
  return(answer_array)
}

## Function to calculate each respondent's score
## as each question is asked
# model - fit model by grm() on test dataset
# answers - answer array returned by adaptive_all()
# n - length of battery
# dataset - test dataset
eval_adaptive <- function(model, answers, n, dataset){
    scores <- errors <- matrix(NA, nrow = nrow(dataset), ncol = n)
    for(i in 1:n){
      current_scores <- factor.scores(object = model, resp.patterns = answers[,,i], method = "EAP")
      scores[ ,i] <- current_scores$score.dat$z1
      errors[ ,i] <- current_scores$score.dat$se.z1
    }
    return(list(scores=scores, errors=errors))
  }

#### Random ####
## Function to calculate one respondent's scores given a random battery
## Returns a 3 dim array, where each slice is a respondent's scores and errors.
# row - used for laply indexing
# model - fit model by grm() on test dataset
# n - length of random battery
# dataset - test dataset
eval_random <- function(row, model, n, dataset){
  random_batt <- sample(colnames(dataset), n)
  scores <- errors <- rep(NA, n)
  
  for(i in 1:n){
    answer_profile <- data.frame(matrix(NA, nrow = 1, ncol = ncol(dataset)))
    colnames(answer_profile) <- colnames(dataset)
    answer_profile[ , random_batt[1:i]] <- dataset[row, random_batt[1:i]]
    current_score <- factor.scores(object = model, resp.patterns = answer_profile, method = "EAP")
    scores[i] <- current_score$score.dat$z1
    errors[i] <- current_score$score.dat$se.z1
  }
  
  return(cbind(scores = scores, errors = errors))
}


#### scaleTester ####
## Function to estimate scores and errors given you know a respondent's
## answers as chosen through a fixed, adaptive, and random survey batteries.
## "True" scores are also simulated.  Function outputs a list of these
## scores, errors, and the GRM models fit to the training and test datasets.
# do_install - should needed packages be installed
# full_data - dataset used for simulation
# num_questions - length of each battery
# num_cores - number of cores to utilize in parallelization
# out_stub - name for output file
# fixed_scale - question names for fixed length battery
scaleTester <- function(do_install = FALSE,
                        full_data,
                        num_questions,
                        num_cores,
                        out_stub,
                        fixed_scale = NULL){
  

  ## Set up
  to_install <- c("ltm", "plyr", "foreach", "doMC")
  if(do_install){
    install.packages(toInstall, repos = "http://cran.us.r-project.org")
  }
  lapply(to_install, library, character.only = TRUE)
  registerDoMC(num_cores)
  
  ## Create training and test datasets
  set.seed(1234)
  trainingSet_size <- floor(nrow(full_data) * 5/6)
  trainingSet_rows <- sample(1:nrow(full_data), trainingSet_size)
  
  trainingSet_data <- full_data[sort(trainingSet_rows), ] 
  testSet_data <- full_data[-trainingSet_rows, ]
  
  ## Fit the models
  ## DOUBLE CHECK -- any more parameters to specify????
  fitModel_training <- grm(trainingSet_data)
  fitModel_test <- grm(testSet_data)
  
  ## Calculate true scores for each person based 
  ## on each full answer profile in test dataset
  true_scores <- factor.scores(object = fitModel_test,
                               resp.patterns = testSet_data,
                               method = "EAP")$score.dat$z1
  
  if(num_questions != length(fixed_scale)){
    stop("num_questions needs to be same length as provided fixed scale")
  }
  
  ## Generating estimates for each battery type ##
  
  ## Fixed
  fixed_results <- eval_fixed(questions = fixed_scale,
                              n = num_questions,
                              model = fitModel_test,
                              dataset = testSet_data)
  fixed_scores <- fixed_results$scores
  fixed_errors <- fixed_results$errors
  

  ## Adaptive
  answer_array <- adaptive_all(model = fitModel_training,
                               n = num_questions,
                               dataset = testSet_data)

  adaptive_results <- eval_adaptive(model = fitModel_test,
                                    answers = answer_array,
                                    n = num_questions,
                                    dataset = testSet_data)
  adaptive_scores <- adaptive_results$scores
  adaptive_errors <- adaptive_results$errors
  

  ## Random
  random_results <- laply(1:nrow(testSet_data),
                          eval_random,
                          dataset = testSet_data,
                          model = fitModel_test,
                          n = num_questions,
                          .parallel = TRUE)
  random_scores <- random_results[,,"scores"]
  random_errors <- random_results[,,"errors"]
  

  sim_output <- list(true_scores = true_scores,
                 adaptive_scores = adaptive_scores,
                 random_scores = random_scores,
                 fixed_scores = fixed_scores,
                 adaptive_errors = adaptive_errors,
                 random_errors = random_errors,
                 fixed_errors = fixed_errors,
                 answer_array = answer_array,
                 fitModel_training = fitModel_training,
                 fitModel_test = fitModel_test)

  save(sim_output, file = out_stub)

  # Mean absolute bias for each battery
  MAB_fixed <- colMeans((sim_output$true_scores - sim_output$fixed_scores)^2
                        + sim_output$fixed_errors^2)
  MAB_dynamic <- colMeans((sim_output$true_scores - sim_output$adaptive_scores)^2
                          + sim_output$adaptive_errors ^2)
  MAB_random <- colMeans((sim_output$true_scores - sim_output$random_scores)^2
                         + sim_output$random_errors^2)
  MAB <- rbind(MAB_fixed, MAB_dynamic, MAB_random)
  rownames(MAB) <- c("Fixed", "Dynamic", "Random")
  
  return(MAB)
}











